Astronomy & Astrophysics manuscript no. AAL7 


February 2, 2008 


(DOI: will be inserted by hand later) 





A method for relaxing the CFL-condition in time explicit schemes 

A. Hujeirat 



Max-Planck-Institut fiir Astronomie, 691 17 Heidelberg, Germany 



cn 
o 
o 

(N 



> 

in 

On ■ 

o : 

m , 
o ■ 

Oh, 

6 ■ 



the date of receipt and acceptance should be inserted later 

Abstract. A method for relaxing the CFL-condition, which limits the time step size in explicit methods in computational fluid 
dynamics, is presented. The method is based on re-formulating explicit methods in matrix form, and considering them as a 
special-Jacobi iteration scheme that converge efficiently if the CFL- number is less than unity. By adopting this formulation, 
one can design various solution methods in arbitrary dimensions that range from explicit to unconditionally stable implicit 
methods in which CFL-number could reach arbitrary large values. In addition, we find that adopting a specially varying time 
stepping scheme accelerates convergence toward steady state solutions and improves the efficiently of the solution procedure. 
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1. Introduction 

The majority of the numerical methods used in astrophysi- 
cal fluid dynamics are based on time-explicit method s (see 
IStone & Norman 1992tEiegler 1998l:lkoide et al. 1999 1) . 
Advancing the solution in time in these methods is based on 
time-extrapolation procedures, that are found to be numerically 
stable if the time-step size is shorter than a critical value, which 
is equivalent to the requirement that the Courant-Friedrich- 
Levy (CFL) number must be smaller than unity. 
This condition, however, limits the range of application and 
severely affects their robustness (see Fig. 1). Using high per- 
formance computers to perform a large number of explicit time 
steps leads to accumulations of round-off' errors that may eas- 
ily distort the propagation of information from the boundaries 
and cause divergence of the solution procedure, especially if 
Neumann type conditions are imposed at the boundaries. 
In this paper we present for the first time a numerical strategy 
toward relaxing the CFL-condition, and therefore enlarging the 
range of application of explicit methods. 



2. Mathematical formulation - scalar case 

In fluid flows, the equation of motion which describes the time- 
evolution of a quantity q in conservation form reads: 



(1) 



where V and / are the velocity field and external forces, respec- 
tively. L represents a first and/or second order spatial operator 
that describe the advection and diffusion of q. 



In the finite space fi, we may replace the time derivative of 



?by: 
~t ~ 



St 



(2) 



where q^ and correspond to the actual value of q at the old 
and new time levels, respectively. 
An explicit formulation of Eq.l reads: 



ot 

whereas the corresponding implicit form is: 



of 

Combining these two approaches together, we obtain: 

§ = e[-LqV + + (1 - 0)[-LqV + f]\ 

ot 



(3) 



(4) 



(5) 



where 0(0 < < 1) is a switch on/off parameter. 

To first order in 6t, we may Taylor-expand q"^^ and /"^' around 



q , I.e., q 

— + 0{LV + g)6q ^ RHS'\ 
St 



q" +Sq + (9((5 ), and obtain: 



(6) 



where g = dfjdq and RMS'" = [/ - LqVf. In the following 
discussion, we term RHS'^ as the time-independent residual. 
Applying Eq. 6 to the whole number of grid points, the follow- 
ing matrix-equation can be obtained: 



(— -I- 0A)5q = RMS'". 
St 



(7) 
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The matrix A contains coefficients such as V/Ax, ri/Ax^ and g 
that correspond to advection, diffusion and to the source terms, 
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Fig. 1. The range of applicability of implicit and explicit methods. The equations describing physical and chemical processes 
whose characteristic time scales are much shorter than the hydrodynamic time scale (thd) should be solved implicitly. For 
example, the characteristic time corresponding to the propagation of gravitational waves (tg), peaks in the radiative energy (tr) 
and chemical processes (tca) occur on much shorter time scales than thd- In most astrophysical problems Alfven wave crossing 
time (tmhd) is shorter than thd- The time scale required for fluid flows (non-relativistic) to relax thermally {tth) is in general 
one order of magnitude larger than thd, depending on the efficiency of the cooling processes. However, the viscous time scale 
(tv,s) can be significantly larger than thd depending on how large the Reynolds number {Re) is. In most astrophysical flows. 
Re is at least 10^. If accretion is considered, the time scale for an envelope to evolve (r^y) is approximately equal to the mass 
of the envelope divided by the accretion rate. Thus, for an envelope of Meii - Mq to evolve from an accretion rate of 
M ~ 10"'*'A1q/F, an accretion time of the order tev ~ IQr' years is required. 
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Fig. 2. The neighboring block matrices in the x and 
y-directions resulting from 5-star staggered grid dis- 
cretization. Entries marked with 'X' denote the elements 
usually used in the impli cit operator splitting approach 
(iHuieirat & Rannacher and '+' are coefficients 

corresponding to the the source terms. The semi-explicit 
method for a scalar equation relies on inverting the diago- 
nal matrix whose entries are marked with X surrounded by 
squares. The generalization of the semi-explicit method to the 
multi-dimensional HD-equations requires inverting the block 
diagonal matrix D^o^- 



respectively. In terms of equation 7, explicit methods are re- 
covered by neglecting the matrix OA. We note that since the 
switch parameter 9 does not appear on the RHS of Eq. 7, ex- 
plicit formulation does not directly depend on 0, but rather on 
the multiplication of OA. Thus, for explicit method to converge, 
it is necessary that the matrix -n OA) is diagonally dominant, 
which implies that OA must be negligibly small. In terms of ma- 
trix algebra, this means that the absolute value of the sum of el- 
ements in each raw of OA must be smaller th an the correspond- 
ing diagonal element 1 16t (jn^^schTliT). In the absence of 
diffusion and external forces (i.e., t] - f - 0), this is equivalent 
to the requirement that at each grid point: 1/dt > |y|/Ax, i.e., 
CFL = 6t\V\/Ax < 1. 

On the other hand, the matrix A can be decomposed as fol- 
lows: A - D + L + U, where D is a matrix that consists of 
the diagonal elements of A. L and U contain respectively the 
sub- and super-diagonal entries of A. Noting that a conserva- 



tive discretization of the advection-difiiision hydrodynamical 
equations (Navier-Stokes equations) gives rise to a D that con- 
tains positive values, we may reconstruct a modified diagonal 
matrix Dmod - I/St + OD. In terms of Eq. 7 we obtain: 



[D^od + 0iL + U)]6q = RHS". 



(8) 



A slightly modified semi-explicit form can be obtained by ne- 
glecting the entries of the matrix 0(L +U).ln this case, a neces- 
sary condition for the iteration procedure to converge is that the 
absolute value of the sum of elements in each raw of 0{L + U) 
must be much smaller than the corresponding diagonal element 
of Dmod in the same raw. In terms of Equation 8, the method is 
said to converge if the entries in each row of D^od fulfill the 
following condition: 



(9) 



roiiowing conaition: 

1 /6t + 0(\V\/Ax + Tj/Ax^ +g)>\\A- D\U 

where ||A - D||co denotes the oo-norm of A - D. This can be 
achieved, however, if the flow is smooth, viscous, and if ap- 
propriate boundary conditions are imposed'. Consequently, the 
inversion process of DmodSq - RHS ° proceeds stably even for 
large CFL-numbers (see Fig. 5). 



3. Generalization 

The set of 2D-hydrodynamical equations in conservative form 
and in Cartesian coordinates may be written in the following 
vector form: 

^+L,,,xf + Ly,yyG=/, (10) 

where F and G are fluxes of q, and Lx.xxj Ly.yy are first and sec- 
ond order transport operators that describe advection-diffusion 
of the vector variables ^ in x and y directions. / corresponds to 
the vector of source functions. 

By analogy with Eq.7, the linearization procedure applied to 
Eq. 10 yields the following matrix form: 



[— -H0(ALx,xx + Biy,yy 
Of 



H)]6q = RHS'', 



(11) 



where A - dF/dq, B - dG/dq and H - df/dq, and which 
are evaluated on the former time level. RHS " = [/ - L^,xxF - 



' Note that diffusion pronounces the inequality in Eq. 9, which gives 
rise to larger CFL-numbers. 
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Adopting a five star staggered grid discretization, it is easy to 
verify that at each grid point Eq. 11 acquires the following 
block matrix equation: 



(12) 



where the underlines (overlines) mark the sub-diagonal (super- 
diagonal) block matrices in the corresponding directions (see 
Fig. 2). D"'^ are the diagonal block matrices resulting from the 
discretization of the operators Lx.xx^', LyjyG and /. 
To outline the directional dependence of the block matrices, we 
re- write Eq. 12 in a more compact form: 



+^^<J^j,k-i- 



(13) 



where Dmod - Sq-^^^l6t+D^ + D^ . The subscripts "j" and "k" de- 
note the grid-numbering in the x and y directions, respectively 
(see Fig. 2). 

Eq. 1 3 gives rise to three different types of solution procedures: 

1 . Classical explicit methods are obviously very special cases 
that are recovered if all the sub- and super-diagonal block 
matrices are neglected, as well as and . The only ma- 
trix to be retained here is (1/5?) x (the identity matrix), i.e., 
the first term on the LHS of Eq. 12. This yields the vector 
equation: 



St 



(14) 



2. Semi-explicit methods are recovered when neglecting the 
sub- and super-diagonal block matrices only, but retaining 
the block diagonal matrices. In this case the matrix equation 
reads: 



(15) 



We note that inverting Dmod is a straightforward procedure, 
either analytically or numerically. 
3. A fully implicit solution procedure requires retaining all the 
block matrices on the LHS of Eq. 13. This yields a global 
matrix that is highly sparse. In thi s case, the "Approximate 
Factorization Method" (-AFM: iBeam & Warming 1973) 
a nd the "Line Gau ss-Seidel Relaxation Method" (-LGS: 
iMacCormack 198^) are considered to be efficient solvers 
for the set of equations in multi-dimensions. 

Fig. 3 shows the time-development of the CFL-number and the 
total residual for 5-different solution procedures for searching 
steady state solutions for Taylor flows between two concen- 
tric spheres. Using spherical geometry, the set of the 2D axi- 
symmetric Navier-Stokes equations are solved. The set consists 
of the three momentum equations, the continuity and the inter- 
nal energy equations. The flow is assumed to be adiabatic. 
In the explicit case, the equations are solved according to 



Eq. 14. For the semi-explicit procedure, we solve the HD- 
equations using the block matrix formulation as described 
in Equation 15. The implicit operator splitting approach is 
based in solving each of the HD-equations implicitly. Here the 
LGS - method is used in the invers ion procedure of each equa- 
tion ( Huieirat & Rannacher 200 H) . Unfortunately, while this 
method has been proven to be robust for modeling compress- 
ible flows with open boundaries, it fails to achieve large CFL- 
numbers in weakly incompressible flows (Fig. 3). This indi- 
cates that pressure gradients in weakly incompressible flows 
do not admit splitting, and therefore they should be included in 
the solution procedure simultaneously on the new time level. 
In the final case, the whole set of HD-equations taking into ac- 
count all pressure terms is solved in a fully implicit manner 
(Fig. 3/bottom). Here we use the AFM for solving the general 
matrix-equation which is locally described by Eq. 13. 
For constructing the time step size in these calcula- 
tions, we have adopted the following description: 6t - 
aoe/ max{RHS j^k), where e = min(AXy, AF^), and where ao 
is a constant of order unity. The maximum and minimum func- 
tions here run over the whole number of grid points. 



3.1. The specially varying time stepping scheme for 
accelerating convergence 

Let [a,b] be the interval on which Eq. 1 is to be solved. We 
may divide [a,b] into N equally spaced finite volume cells: 
Axi - (b - a) IN, i = 1, A^. To follow the time-evolution of 
q using a classical explicit method, the time step size must ful- 
fill the CFL-condition, which requires 5t to be smaller than the 
critical value: = min{Axi/(y + Vs)i}. 
If [a,b] is divided into N highly stretched finite volume cells, for 
example Axi < Ax2... < Axn, then the CFL-condition restrict 
the time step size to be smaller than St"" - min,{(A.Xi/(y-i-ys)i}, 
which is much smaller than 5f". Thus, applying a conditionally 
stable method to model flows using a highly non-uniform dis- 
tributed mesh has the disadvantage that the time evolution of 
the variables in the whole domain are artificially and severely 
affected by the flow behaviour on the finest cells. 
Moreover, advancing the variables in time may stagnate if the 
flow is strongly or nearly incompressible. In this case, Vs » V, 
which implies that the time step size allowed by the CFL- 
condition approaches zero. 

However, we may associate still a time step size with each grid 
point, e.g., 5f™ = Axi/(y + Vs), and follow the time evolu- 
tion of each variable qi independently. Interactions between 
variables enter the solution procedure through the evaluation 
of the spatial operators on the former time level. This method, 
which is occasionally called the "Residual Smoothing Method" 
proved to be efficient at providing quasi-stationary solutions 
within a reasonable number of iteration, when compar ed to 
normal explicit methods. (Fig. 3, also see 'Enand er 19971) . The 
main disadvantage of this method is its inability to provide 
physically meaningful time scales for features that possess 
quasi-stationary behaviour. Here we suggest to use the obtained 
quasi-stationary solutions as initial configuration and re-start 
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the calculations using a uniform and physically well defined 
time step size. 

4. Summary 

In this letter we have presented a strategy for relaxing the CFL- 
condition which enlarges the range of application of explicit 
methods and improve their robustness. The method is based on 
re-formulating explicit methods in matrix-form, which can be 
then gradually modified up to a fully implicit scheme. The ma- 
trix corresponding to the semi-explicit scheme presented here 
is a block diagonal matrix that can be easily inverted, either 
analytically or numerically. Unlike normal explicit methods, 
in which the inclusion of diffusion limits further the time step 
size, in the semi-explicit formulation presented here diffusion 
pronounces the diagonal dominance and enhances the stabil- 
ity of the inversion procedure, irrespective of the dimensional- 
ity of the problem. We note that the CFL-numbers achieved in 
the present modeling of Taylor flows are, indeed, larger than 
unity, but they are not impressively large as we have predicted 
theoretically. We may attribute this inconsistency to three dif- 
ferent effects: 1) The flow considered here is weakly incom- 
pressible. This means that the acoustic perturbations have the 
largest propagation speeds, which requires that all pressure ef- 
fects should be included in the solution procedure simultane- 
ously on the new time level. 2) The conditions imposed on the 
boundaries are non-absorbing, and do not permit advection of 
errors into regions exterior to the domain of calculations. 3) 
The method requires probably additional improvements in or- 
der to achieve large CFL-numbers. This could be done, for ex- 
ample, within the context of the "defect-correction' iteration 
procedure, in which the block diagonal matrix D^od is em- 
ployed as a pre-conditioner. 

Finally, we have shown that the residual smoothing ap- 
proach improves the convergence of explicit methods. 
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Fig. 3. The development of the CFL-number (left axis) and the 
total residual (right axis) versus covered-time in normalized 
units of five diff'erent numerical methods (from to top to 
bottom: normal explicit, residua l smoothing, semi-exp licit, 
implicit operator splitting (.Huieirat & Rannacher 20011) and 
the fully implicit method). While the eff'ective time covered 
in each run of these different methods is similar, the actual 
number of iteration is substantially diff'erent. The numerical 
problem here is to search stationary solutions for Taylor flow 
between two concentric spheres. The inner sphere has a radius 
r,„ - 1 and rotates with angular velocity = 5 , whereas 
the outer sphere is non-rotating and its radius is taken to be 
rout = 1.3 . We use the viscosity coefficient v- 10"^. The initial 
density and temperature are taken to be p(r, 6*, f = 0) = 1 , and 
T{r,6,t = 0) = 10', respectively. The computational domain 
is [1,1.3] X [0,7r/2] and consists of 30x50 non-uniformly 
distributed tensor-product mesh. 
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Fig. 4. Steady state solutions of the Taylor flow between two 
concentric spheres. Here the velocity field and the angular ve- 
locity (large-to-low values correspond to blue-to-red colors) are 
shown. The capability of the methods to capture the formation 
of rotationally-driven multiple vortices near-equatorial region 
is obvious. 
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Fig. 5. The semi-explicit method applied to a scalar problem. 
The Figure shows the development of the CFL-number (left 
axis) and the residual (right axis) of the angular momentum 
equation in two dimensions versus the iteration number As ini- 
tial conditions we use the steady distributions of the physical 
variables that have been obtained from the simulations of the 
Taylor problem (see Fig. 4). This includes the velocity field, 
density, temperature and 77. For the angular velocity we use 
Q = as initial condition. 



